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Abstract 

This article introduces the Parabolic Variance (PVAR), a wavelet variance similar to the Allan variance, based on the Linear 
Regression (LR) of phase data. The companion articl^arXiv: 1506.05009 [physics.ins-det] details the Q frequency counter, which 
implements the LR estimate. 

The PVAR combines the advantages of AVAR and MVAR. PVAR is good for long-term analysis because the wavelet spans 
over 2t, the same of the AVAR wavelet; and good for short-term analysis because the response to white and flicker PM is l/r^ 
and l/r^, same as the MVAR. 

After setting the theoretical framework, we study the degrees of freedom and the confldence interval for the most common 
noise types. Then, we focus on the detection of a weak noise process at the transition - or corner - where a faster process 
rolls off. This new perspective raises the question of which variance detects the weak process with the shortest data record. Our 
simulations show that PVAR is a fortunate tradeoff. PVAR is superior to MVAR in all cases, exhibits the best ability to divide 
between fast noise phenomena (up to flicker FM), and is almost as good as AVAR for the detection of random walk and drift. 


1. Introduction 

The Allan variance (AVAR) was the first of the wavelet-like variances used for the characterization of oscillators and 
frequency standards (Tj. After 50 years of research, AVAR is still unsurpassed at rendering the largest r for a given time series 
of experimental data. This feature is highly desired for monitoring the frequency standards used for timekeeping. 

Unfortunately, AVAR is not a good choice in the region of fast noise processes. In fact, the AVAR response to white and 
dicker PM noise is nearly the same, 1/r^. For short-term analysis, other wavelet variances are preferred, chiefiy the modified 
Allan variance (MVAR) H) 0 The MVAR response is 1/r^ and 1/r^ for white and dicker PM, respectively. However, 
MVAR is poor for slow phenomena because the wavelet spans over 3r instead of 2r. Thus, for a data record of duration T, 
the absolute maximum r is T/3 instead of T/2. 

Speaking of ‘wavelet-like’ variances, we review the fundamentals. A wavelet is a shock with energy equal to one and 
average equal to zero, whose energy is well confined in a time interval (see for example p. 2]) called ‘support’ in proper 
mathematical terms. In formula, dt = 1, dt = 0, and dt = ^ ~ with small e > 0. It makes 

sense to re-normalize the wavelet as ^ dt = 1, so that it is suitable to power-type signals (finite power) instead of 

energy-type signals (energy finite). By obvious analogy, we use the terms ‘power-type wavelet’ and ‘energy-type wavelet’. 
These two normalizations often go together in spectral analysis and telecom (see the classical books ||^, 0 )- For historical 
reasons, in clock analysis we add a trivial coefficient that sets the response to a linear drift Dy to |Dy , the same for all the 
variances. 

High resolution in the presence of white and dicker phase noise is mandatory for the measurement of short-term fiuctuations 
(/iS to s), and useful for medium-term fiuctuations (up to days). This is the case of optics and of the generation of pure 
microwaves from optics. The same features are of paramount importance for radars, VLBI, geodesy, space communications, 
etc. As a fringe benefit, extending the time-domain measurements to lower r is useful to check on the consistency between 
variances and phase noise spectra. MVAR is suitable to the analysis of fast fiuctuations, at a moderate cost in terms of computing 
power. Frequency counters specialized for MVAR are available as a niche product, chiefiy intended for research labs 0 . 

A sampling rate of 1/r is sufficient for the measurement of AVAR, while a rate of I/tq = m/r is needed for MVAR, 
where the rejection of white phase noise is proportional to m. MVAR is based on the simple averaging of m fully-overlapped 
(spaced by the sampling step tq) frequency data, before evaluating 

The linear regression provides the lowest-energy (or lowest-power) fit of a data set, which is considered in most cases as 
the optimal approximation, at least for white noise. For our purposes, the least-square fit finds an obvious application in the 
estimation of frequency from a time series of phase data, and opens the way to improvements in fiuctuation analysis. Besides, 
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new digital hardware — like Field-Programmable Gate Arrays (FPGAs) and Systems on Chip (SoCs) — provides bandwidth 
and computing power at an acceptable complexity, and makes possible least-square fitting in real-time. 

We apply least-square estimation of frequency to fast time stamping. The simplest estimator in this family is the linear 
regression (LR) on phase data. The LR can be interpreted as a weight function applied to the measured frequency fluctuations. 
The shape of such weight function is parabolic. The corresponding instrument is called 'Q counter,’ described in the companion 
article |[^. The name ft comes from the graphical analogy of the parabola with the Greek letter, in the continuity of the If 
and A counters p0| , (TT). The Q estimator is similar to the A estimator, but exhibits higher rejection of the instrument noise, 
chiefly of white phase noise. This is important in the measurement of fast phenomena, where the cutoff frequency fn is 
necessarily high, and the white phase noise is integrated over the wide analog bandwidth that follows. 

In the same way as the Ft and A estimators yield the AVAR and the MVAR, respectively, we define a variance based on 
the Q estimator. Like in the AVAR and MVAR, the weight functions are similar to wavelets, but for the trivial difference that 
they are normalized for power-type signals. A similar use of the LR was proposed independently by Benkler et al. (H at the 
IFCS, where we gave our first presentation on the O counter and on our standpoint about the PVAR. In a private discussion, 
we agreed on the name PVAR (Parabolic VARiance) for this variance, superseding earlier terms (H). 

We stress that the wavelet variances are mathematical tools to describe the frequency stability of an oscillator (or the 
fluctuation of any physical quantity). Albeit they have similar properties, none of them should be taken as “the stability” of 
an oscillator. For the same reason, MVAR and PVAR should not be mistaken as ‘estimators’ of the AVAR. To this extent, the 
only privilege of AVAR is the emphasis it is given in standard documents (E)- 

After setting the theoretical framework of the PVAR, we provide the response to noise described by the usual polynomial 
spectrum. Then we calculate the degrees of freedom and confidence intervals, checking on the analytical results against extensive 
simulations. Finally, we compare the performance of AVAR, MVAR and PVAR for the detection of noise types, using the 
value of r where cr^(r) changes law as an indicator. In most practical cases PVAR turns out to be the fastest, to the extent 
that it enables such detection with the shortest data record. 

IT Statement oe the Problem 

The clock signal is usually written as 

v{t) = Vq sin[27rz/of + ^{t)] 

where Vq is the amplitude, z/q is the nominal frequency, and ip{t) is the random phase fluctuation. Notice that ip{t) is allowed 
to exceed ±7r. Alternatively, randomness is ascribed to the frequency fluctuation (Az/)(t) = 27rip{t). 

We introduce the normalized quantities 

x(t) = t x{t) 

y{t) = 1 + y(i) 

where x(t) = (/2(t)/27rz^o, and y(t) = x(t). The quantity x(t) is the clock readout, which is equal to the time t plus the random 
fluctuation x(f). Accordingly, the clock signal reads 

v{t) = Vq sin[27rz/ox(0] 

= Vb sin[27rz/o^ + 27ruQx{t)] 

For the layman, x is the time displayed by a watch, t is the ‘exact’ time from a radio broadcast, and x the watch error. The error 
X is positive (negative) when the watch leads (lags). Similarly, y is the normalized frequency of the watch’s internal quartz, 
and y its fractional error. For example, if y = -blO ppm (constant), the watch leads uniformly by 1.15 s/day. For the scientist, 
x{t) is the random time fluctuation, often referred to as ‘phase time’ (fluctuation), and y(t) is the random fractional-frequency 
fluctuation. The quantities x(t) and y(t) match exactly x{t) and y{t) used in the general literature of time and frequency 

GD GD 

The main point of this article is explained in Fig. We use the linear regression of phase data to get a sequence {y} of 
data averaged on contiguous time intervals of duration r, and in turn the sequence {y} of fractional-frequency fluctuation data. 
Two contiguous elements of {y} and {y} are shown in Fig. from which we get one value of ^(y 2 — yi)^ for the estimation 
of the variance. 

Most of the concepts below are expressed in both the continuous and the discrete settings with common notations without risk 
of confusion. For example, the same expression x = t-\-x maps into =tk-\-Xk in the discrete case, and into x{t) = t-\-x{t) 
in the continuous case. The notations (.), (.,.) and || . || represent the average, the scalar product and the norm. They are 
defined as (x) = ^ (^^v) = n 11^11 = (x, x)^/^ where n is the number of terms of the sum in the discrete 

case, and as (x) = ^ f x(t) dt^ (x, y) = ^ f x(t)y(t) dt, ||x|| = (x, x)^/^ where T is the length of the interval of integration in 
the continuous case. The span of the sum and the integral will be made precise in each case of application. The mathematical 
expectation and the variance of random variables are denoted by IE{ . } and V{ . }. 
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Fig. 1. Principle of two-sample linear regression measurement, and notation. 


The linear regression problem consists in searching the optimum value y of the slope 77 (dummy variable) that minimizes 
the norm of the error x — xq — rjt, i.e., y = argmin^ ||x — xq — r]t\\‘^. Since we are not interessed in xq, which only reflects 
choice of the origin of x, the solution is the random variable 

^ (x- (x) 

^ \\t-{t)\\^ ■ 

We recall some useful properties of y as an estimator of the slope of x. For the sake of simplicity, with no loss of generality, 
we refer to a time sequence is centered at zero, i.e., (t) = 0 . 

1) The estimator y can be simplifled as 

y=M 


2) If the component (or the values x(t)) are independent, the estimator variance is 

v{y} = 


The assumption of independent continuous random process is rather usual in theoretical works. However this is done to 
simplify some proofs, the results can be used in their discrete form. 

3) Sampling uniformly at the interval tq, the discrete time is tk = {k ^ |)ro for k G {—p, m = 2p and r = mro. 

For large m, we get 

y«H-and V{y} «- 

4) With a signal that is continuous over a symmetric time interval (—^, p, we get 

12{x,t) 


1 + 


( 1 ) 


The continuous form of the estimator y can be expressed as a weighted average of x or y. For this purpose, it is useful to take 
y as a time dependent function deflned over t G ( 0 , r) 

12 

^ “3 / sx{t-r/2 + s) ds 

^ J-t/2 


12 r 

= (s + r/ 2 ) x(t + s) ds 

f {r/2 — s) x{t — s) ds. 

Jo 


12 


( 2 ) 


III. Time Domain Representation 

A. Generic Wavelet Variance 

Let us denote with T the duration of the data run, with tq the sampling interval, with N the number of samples, and with 
n the ratio T/r. Thus, T = Ntq, and N = mn. We consider the series {yi}i=i,..,n of frequency deviation estimates. In this 
section we denote with cr^(r) a generic wavelet variance, either AVAR, MVAR, PVAR, etc. 
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In the case of uncorrelated frequency fluctuations (white FM), an unbiased estimator of the variance V {y} is 

n ^ n ^ 




i=i 


so 


V{y}=E{52_i}. 

After Allan Q, we replace the estimator with a two-sample variance by setting n = 2. Then, the variance V {y} = E {5'^} 
is 

= 2® {(yi “y2)^}, 


and its estimator averaged over the n — 1 terms 


\t) = i ((y* -Yi+i)^) . 


Notice that two-sample variance is generally written as (r), and that we drop the subscript y. 
Following the Lesage-Audoin approach (H) we deflne the point variance estimates 


and the estimated variance 






(3) 


(4) 


(5) 


( 6 ) 


i=l 


The relationship between the ai and the N individual x{kTo) measures depends on the type of counter (Ft, A, O). 


B. Continuous-Time Formulation of PVAR 

In the case of continuous time, the difference between contiguous measures is 


N N 12 

y{t + r)-y{t) = 


12 


12 


I {l-s)<t + r-s)ds- (Z-s)xit-s) 


ds 


— J {\t - s\ - f) y<{s) ds. 


Accordingly, the two-sample variance 0 is written as 






and notice the subscript P for PVAR. Such variance is independent of t, and it can be expressed as the running average 

2 \ 


where 


is the even weight function, and 


crp(r) = E<| ^ J x{s) Wy^{s — t) 
6^2/^,,, r\ 

2) A(-r,r)(^) 


(7) 


UJx(t) = 


X( — t,t) (l') 


1 te{-r,r) 
0 elsewhere 


is the indicator function (or characteristic function). 

From (0, we see that PVAR can also be written as a convolution product 

2> 


crp(r)=E|^y x(s) f)x(i - s) | = e| ^x(t) * f)x(i) 


where ^{t) is the convolution kernel which applies to x(f). The kernel ()(t) is related to the weight function w{t) by the general 
property that ()(t) = w{—t). However, since Wx{t), is even function, it holds that ^x{t) = Wx{t). 
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Fig. 2. Convolution kernels of PVAR from phase data (above) or frequency deviations (below) for r = Stq. 


Similarly, the estimator 0 is written 


as 


where 


= E 




3V2i ,, , 


1 2 


dt 


) X( —T,r) (0 


is the convolution kernel which applies to y(t). 

Thanks to the fact that y{t) = x(f), crp{r) can also be expressed as the running average 

2 > 


aj,{T) =Ei^^J^y{s)Wy{s-t) ds^ |, 


( 8 ) 

(9) 


where 

Wy{t) = -r) X(-T,r){t) 

is the weight function. Since Wy{t) is odd function, it hold that ()y(t) = —Wy(t). Moreover, the parabolic shape of the PVAR 
wavelet comes from the t • |t| factor in Wy{t) and f)y(t). 

For the purpose of operation with the Fourier transform, it is convenient to restate these expression in terms of filter or 
convolution 


The weight functions Wy^(t) and Wy{t), and also the kernels ^x{t) and i)y(t), match the definition of power-type wavelet given 
in the introduction. As a consequence of the property y(t) = x(t), it holds that ^x{t) = Figure^ shows the convolution 

kernels associated to PVAR. 

It is worth pointing out that our formulation is is general, as it applies to AVAR, MVAR, PVAR, and to other similar 
variances as well. Of course, the wavelet depends on the counter (Fig [^. 
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Fig. 3. The wavelets associated to AVAR, MVAR, and PVAR. 


C Practical Evaluation of PVAR 

Denoting the discrete time with = x(iro), the estimate of the two-sample variance is Gz) 

for AVAR, with M = N — 2m, and 

m—1 


ai = 


V2 


^ ^ ( ^i+k A' ^^i+m+k Xi+2m+/cl 
1 mr ^ V / 


/c=0 


( 11 ) 


( 12 ) 


for MVAR, with M = N — 3m + 1. 

Now we calculate for PVAR. First, the discrete form of y can be obtained from ^ by replacing the time integral with 
a sum with a time increment equal to tq. Accordingly, r is replaced with mro, s with kro, t with iro, and x{t) with x^ 


Similarly, 


and consequently 


Yi = 


12 

m^TQ 


m—l 


E 


/ (m - 1 )to 

V 2 


- kro 


^i—k'^O 


12 


m—l 

E 


/c=0 


m ■ 



^i—k 


Yi+i 


-I o 3 / I 

12 \ / 777 / — 1 


/c=0 


^ ) ^i+m—A 


Yi - Yi+i = 

m^T 


12 ^ 777 / — 1 ^ ^ ^ 

> , I -^- k I [Xi-k - ^i+m-k) • 


/c=0 
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Second, we recall that is defined for i > 0. Hence, we have to shift the origin by m — 1, so that also is defined with 

i = 0 


12 

Yi - Yi+i = 


m-l , ^ ^ X 

^ ^ f ^ ^ j — k Xi+2m—1 —/c) • 

/c=0 ^ 


Third, since the coefficient (m — l)/2 — /c is symmetrical for k and for m — 1 — k, we interchange i m — 1 — k with 
i m — 1 — {m — 1 — k) = i k, and i + 2m — 1 — k with i + 2m — 1 — {m — 1 — k) = i m k 


12 

Yi - Yi+i = 

m^T 


m-l y ^ ^ X 

; 2 (Xz + Zc ~ + + . 


Finally, it comes 


6a/2 


m-l ^ X 

^ f ^ - kj (xi+k - Xi+m+/c) for m > 2 

k=0 ^ ' 


OLi = 

M = — 2m + 2. 


(13) 


For consistency with AVAR and MVAR, we require crp(ro) = cr^(ro) = cr|^(ro), i.e. all variances are equal at sampling time 
tq. Since •EH* gives = 0 for m = 1, we redefine 


a. 


V2to 
M = N-2 


— ^ “h 2 x2-|-1 X2-|-2^ 


for m = 1 


Having N samples {x^} taken at the interval r, the estimate d-p{r) can be computed using and (13) as 


aUr) = 


72 




M-l 

E 

i=0 


'm-l y ^ _ 2 \ ' ^ 

./c=0 ^ ^ 


(14) 


(15) 


D. Time-Domain Response 
From it comes 

cr|,(r) = E{,7 |,(t)} 


72 


M-l 


E« 


777,4^2 I ^ 


- E 


i =0 


'm—1 


( O ^ + ^ + + 

/c=0 




m-l m-l / 1 \ / 1 

72 ( m — 1 \ ( m — 1 

k][ —-- I 


EE 


k=0 1=0 

Rx{{k - l)ro) - Rx{{k -m- l)ro) - Rx{{m ^k- l)ro) + Rx{{k - l)ro) 


(16) 


where Rx{0) = E {x(t)x(t + 6>)} is the autocorrelation function of x(t). The autocorrelation function is detailed in Section 


V-D Whereas Rx{r) depends on /l and fn, these parameters cancel in the derivation of PVAR. 


IV. Frequency Domain Representation 


A. Transfer Function 

The transfer function Hp{f) of PVAR is the Fourier transform of the kernel ()y(f). The square of its modulus is given by 


\Hp{f)f = 


9 [2sin^(7rr/) — ttt/ sin(27rr/)] ^ 
2(7rr/)6 


(17) 


Figure 1^ shows |i7p(/)|^, together with the transfer function of AVAR and MVAR. All are bandpass functions with approxi¬ 
mately one octave bandwidth. However, PVAR exhibits significantly smaller side lobes because the weight function is smoother. 
This is well known with the taper (window) functions used in the digital Fourier transform p~8| . 





























Fig. 4. PVAR transfer function, compared to AVAR and MVAR, for integration time is r = 1 s and sampling interval tq = t I^ = 250 ms. 


This can be proved as follows. The transfer function is obtained after Fourier transformation, using the property that ()y 
odd function 


Hy{f) = [ 

JR 


\/2t 

The primitive is calculated by parts integration 

3i 


le 


— 2i'KTf 




Then, 


Hp{f) 




1 — cos (27rr/) — ttt/ sin (27rr/) 


Finally, using 1 — cos (27rr/) = 2sin^ (ttt/), we get 


and 


Hp{f) = 


3i 


9 


2 sin^ (ttt/) — ttt/ sin {^nrf) 


27r6^6j6 _ 


1 — cos (27rr/) — ttt/ sin (27rr/) 


1 2 


B. Convergence Properties 
For small /, it holds that 

sin(7rr/) « ttt/ - + 0{P) 

sin(27rr/) « 27rr/ - + 0{f) 

SO 

2 sin^ (ttt/) — ttt/ sin (27rr/) « 2 ^ttt/—— ttt/ ^27rr/ 

« + 12)+fOU) 

then, at low frequency, 

Hp{f) « 'Pli'KTf. 

We conclude that 

|JTp(/)p ^ 27r^r^/^ at low frequency, 
thus PVAR converges for 1//^ FM noise. Similarly 

\Hp{f)\^ oc (ttt/)”^ at high frequency, 
therefore PVAR converges for FM noise. 


4 

3 


TT r / 































TABLE I 

Response oe AVAR, MVAR and PVAR to the common noise types, and to driet. 


Noise 

type 

SyU) 

AVAR (7\{t) 

MVAR alfir) 

PVAR 


White 

PM 


3h2 

3h2 

3h2 

4 





Flicker 

PM 

hi/1 

[1.038 + 3 ln(7rT/To)] hi 

47r2^2 

[241n(2) -91n(3)]hi 

3 [ln(16) - 1] hi 

3.2 

White 

FM 

ho/° 

ho 

2r 

ho 

4t 

3ho 

5t 

2.4 

Flicker 

FM 


21n(2)h_i 

[271n(3) - 321n(2)] h_i 

2[7-ln(16)]h_i 

1.8 

8 

5 

Random 

walk FM 

h-2/-2 

27r^h 2T 

ll7r^h 2T 

267r^h_2T 

1.4 

3 

20 

35 

Drift 

y(t) = Dyt 




1 

The lowpass cutoff frequency fn, needed for AVAR, is set to 1/2to (Nyquist frequency) 



Fig. 5. Response of PVAR to the polynomial-law noise types, and to linear drift. 

C. Calculation of PVAR from Spectral Data 

Given the Power Spectral Density (PSD) S^{f), PVAR evaluated as 

nOO 

a%{r)= \Hp{ffSy{f)df. 
Jo 
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Replacing Sy{f) with the components of the polynomial law, from h-2f~^ (random walk FM) to h2p (white PM), we get the 
response of shown on Table |Ij together with AVAR and MVAR. Figure shows the response of PVAR to the polynomial-law 
noise types as a function of the integration time r. 

For comparison, crp(r) can also be calculated in the time domain using ( p^ , and also with Monte-Carlo simulations (see 
Section [V-E| ). Time domain, frequency domain and Monte-Carlo simulations give fully consistent results. 


V. Degrees oe Freedom and Confidence Interval 


A. Equivalent Degrees of Freedom (EDF) 

We consider the estimates of a generic variance cr^(r), assumed -distributed, k G The EDF a depend on the 
integration time r, and of course on the noise type. The mean and variance (the variance of the variance) are 


= kE{xl} = kv 
V{6-2(r)} = PV{xl] = 2k^v. 
Accordingly, the degrees of freedom u are given by 

^ 2E{a^{T)}^ 

^ V{<T2(r)} • 


(19) 


Thus, the knowledge of a enables to define a confidence interval around E {(3'^(t)} with given confidence p. For applying this 
result to PVAR, we have then to calculate the variance of PVAR. 


B. Variance of PVAR 

The variance of the estimate (J^(r) is given by 

V{6-2(r)} =e|[6-2(t)-E{CT^( r)}]^| 


Expanding ( [20| yields 

i =0 i =0 

The Isserlis’s theorem (2D states that, for centered and jointly Gaussian random variables ^ and w 

-E{m;2} =2[E{^^<;}]^ 

Assuming that X is a Gaussian process and that ai, aj are two centered jointly Gaussian random variables, it comes 

M-lM-l 2 

i=0 j=0 

The derivation of E{aiaj} is given in the next Section. 


= E 


M-l ( M-1 

i=0 V i=0 


( 20 ) 


( 21 ) 


C. Equivalent Degrees of Freedom of PVAR 
From ( p^ , we can calculate 

¥.{aiaj} = 


72 




E« 


m—1 


( O ^ Xi+m+/c) 

lk=0 


'm-l / _ -1 \ 

E - y ~ ^j+m+l) 


which expands as 


E{aiaj} = 




L 1=0 


m-l m-l X 1 \ / 1 

72 ( m — 1 ^\ ( m — 1 


EE 

k=0 1=0 


|2i?x[(i + k-j- 1 )to] 


-i?x[(« + k- j-m- 1 )to\ - i?x[(i + m + k- j - Ot-q]!- 


(22) 
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S.if) 


TABLE II 

Autocorrelation function of the phase-time fluctuation. 


R.jO) 


R,{t) (for T / 0) 


ko 

k-i/-' 

k-2/-2 

k-3/-" 


k-4/ 


k-i 


ko/if 


- +\n{fH/fL) 


k-2 


k-3 


1 1 
./l fn. 


1 1 
72 “ 7^ 

■/L J H 


k-4 

3 


1 1 
73 “ 7^ 

■/L J H 


k-1 


0 


cos(27r/Lr) - 1 + 27t f lt sui{27T f lt) 


(27r/LT)2 

+ Ci(27rT/if) - Ci(27rT/L) 


k-2 


cos(27rfLT) cos{27rfHT) 


II fn 

+27rT [si(27r/z,T) - Si(27r/ifT)j | 


k-3 


COs(27vfLr) COs{27TfHT) 


vl vl 

+27r^T^ |^Ci(27r/LT) - Ci(27r/ifT)j 
sin(27r/ifT) sin(27r/Lr) 


+7rT 


l^sir 


fn 


fL 


]} 


k_4 


J (27r^/^T^ - 1) cos(27r/jf t) + TTfar sin(27r/jf t) 

1 

{27r^flT^ - 1) cos{27rf lt) + tt/^t sin(27r/LT) 
^47 i^t _ Si(27r/LT)] 


fL and fn are the highpass and lowpass cutoff frequencies which set the process bandwidth 
Ci(x) and Si(x) are the Cosine and Sine Integral functions 


Thanks to (21) and (22), we can calculate the variance of PVAR from the autocorrelation function Rx{r). For example, in the 
case of white PM noise we find 


and 


V {o-|,(t)} 


9hi 


m 

23t7 - 12 
M 



- 175 


M2 


35 

23m/M - 12(m/M)2 - 175m/M2' 


(23) 

(24) 


D. Numerical Evaluation of the EOF 

The EOF can be evaluated by substituting ( [22l > into ( |^ , and then into ( p^ . In turn, thanks to the Wiener Khinchin 
theorem, stationary ergodic processes states that Rx{r) can be obtained as the inverse Fourier transform of the Power Spectral 
Density (PSD). Since the PSD is real and even p^ , p3| , we get 

E~\-00 

Rx{t)= Sx{f) cos{27rrf) df. (25) 

Jo 

Replacing Sx{f) with the polynomial law from white PM to random walk FM (f~^ PM), we get the results shown in Table [n| 
The derivation is rather mechanical, and done by a symbolic algebra application (Wolfram Mathematica). For numerical 
evaluation — unless the reader understanding the computer code in depth — we recommend the approximations lim^o Ci(x) = 
C + ln(x) — x^/96, where C ^ 0.5772 is the Euler-Mascheroni constant, lima^^oo Ci(x) = 0, lima^^o Si(x) = x — x^/9, and 
liiHa^^Qo Si((r) — 'KI 2 . 

As an example, we take a data record of = 2048, tq = 1 s, high cut-off frequency fn = ^ (equal to the Nyquist 
frequency), low cut-off frequency /l = 256 Atq physical meaning of /l) and r G {tq, 2ro, 4ro,..., 1024ro}. 

Figureshows the EDF for the common noise typesl^ooming in (Fig. [fright), we see that the plots do not overlap. 

E. Monte-Carlo simulations 

Another way to assess the EDF is by simulated time series. We generated 10000 sequences of = 2048 samples for each 
type of noise using the 'RruiteuC' noise simulator p4| , which is based on filtered white noise. This code is a part of the 
SigmaTheta software package, available on the URL given by p4| . It has been validated by more than 20 years of intensive 
use at the Observatory of Besancon. Again, the EDF are calculated using 
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Fig. 6. Numerical computation of the PVAR EDF for the common noise types. The right hand plot is a crop of the left one. 



Fig. 7. Comparison of the EDF calculated analytically by numerical computation (Sec. [V-D^ , and assessed by Monte-Carlo simulation (Sec. [V-E^ . 

In the end, we compared three methods, the autocorrelation function Rx{r) with and ( [^ , the Monte-Carlo simulation 
with hruiteur code, and the analytical solution ( [24| ), the latter only with white noise. Figure [7] compares the EDF obtained with 
these three methods. The results match well, with a discrepancy of a few percent affecting only the first two points (r < 2ro). 
The reason is that, with such a small t/tq ratio, the Wy weight function is a poor approximation of the parabola of the Q 
counter (see j^). 

Table and Fig. compare the EDF of PVAR to AVAR and MVAR. MVAR is limited to r = 682 because the wavelet 
support (span) is 3r instead of 2r. 


VI. Detection of Noise Processes 

Running an experiment, we accumulate a progressively larger number N of samples x/.. As N gets larger, we fill up the 
(jy (r) plot adding new points at larger r. Besides, at smaller r the error bars shrink because the number of degrees of freedom 
increases. Looking at the log-log plot, we find the fast processes on the left and the slow processes on the right. This is due to 
the nearly-polynomial law of Table Having said that, we address the question of which variance is the most efficient tool 
at detecting a slower process ‘SP’ in the presence of a faster process ‘FP’ as illustrated in Fig.|^ The criterion we choose is 
the lowest level of the SP that can be detected 

• with a probability of 97.5% (i.e., two sigma upper bound) 

• in the presence of the faster process FP of given level, 

• using a data record of given length N. 

Our question about the most efficient tool relates to relevant practical cases detailed underneath. 

Our comparison is based on a simulation with N = 2048 samples uniformly spaced by tq = 1 s. So, the lowest r is equal 
to 1 s, and the largest r is equal to Nto/2 = 1024 s for AVAR and PVAR, and to Ntq/S = 682 s for MVAR. 

For fair comparison, we re-normalize the variances for the same response to the SP process. For example, the response 
to white FM noise Sy{f) = ho is ho/2r for the AVAR, ho/4r for the MVAR, and 3ho/5r for the PVAR. Accordingly, a 












EDF EDF EDF 


13 


TABLE III 

Comparison of the EDF of AVAR, MVAR and PVAR for the common noise types. 


r/ro 

1 

2 

4 

8 

16 

32 

64 

128 

256 

512 

682 

1024 

White PM (/+^ EM) 











AVAR 

892 

1060 

1020 

1010 

955 

953 

922 

896 

811 

652 


0.981 

MVAR 

891 

970 

685 

355 

173 

82.5 

38.9 

17.3 

7.48 

2.88 

1.02 


PVAR 

892 

1150 

824 

419 

202 

99.1 

46.9 

22.0 

10.0 

4.13 


1.03 

Flicker PM (/+^ EM) 











AVAR 

1090 

1140 

984 

728 

523 

340 

209 

127 

69.5 

33.8 


0.930 

MVAR 

1090 

1020 

544 

258 

126 

62.1 

29.3 

13.9 

5.73 

2.09 

1.04 


PVAR 

1090 

1300 

701 

329 

165 

79.4 

38.2 

18.4 

8.42 

3.36 


1.05 

White EM (/^ EM) 











AVAR 

1380 

1200 

716 

372 

186 

91.7 

45.3 

21.8 

10.2 

4.07 


1.01 

MVAR 

1380 

1060 

505 

247 

119 

58.4 

28.6 

13.2 

5.71 

1.87 

1.04 


PVAR 

1380 

1390 

680 

319 

157 

76.7 

37.5 

18.2 

8.43 

3.32 


1.01 

Flicker EM 

EM) 











AVAR 

1780 

1200 

595 

299 

150 

72.8 

36.1 

17.1 

7.58 

3.05 


1.02 

MVAR 

1780 

1030 

484 

241 

120 

57.9 

28.5 

12.9 

5.32 

1.58 

1.02 


PVAR 

1780 

1470 

648 

319 

159 

77.8 

38.2 

18.2 

8.01 

3.16 


1.02 

Random 

walk EM (/-■"= 

EM) 










AVAR 

1990 

1020 

480 

238 

117 

57.9 

28.1 

13.3 

5.93 

2.29 


1.01 

MVAR 

1990 

861 

398 

197 

96.5 

47.1 

22.6 

10.3 

4.26 

1.31 

1.02 


PVAR 

1990 

1290 

548 

266 

131 

64.3 

31.2 

14.8 

6.53 

2.49 


1.02 


1e+03 


1e+02 


1e+01 


1e+00 


1e+03 


1e+02 


1e+01 


1e+00 


1e+03 


1e+02 


1e+01 


1e+00 

1e+00 1e+01 1e+02 1e+03 

Integration time t 




1e+00 1e+01 1e+02 1e+03 

Integration time x 


Fig. 8. Comparison of the Equivalent Degrees of Freedom of AVAR, MVAR and PVAR for the different types of noise. All noise sequences were simulated 
with a unity coefficient noise, 2048 samples and a sampling frequency of 1 Hz. 
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Fig. 9. Concept of noise process detection. The process SP is barely visible in A (pe50% probability), detected with a probability of 97.5% in C (threshold 
of nearly certain detection), and detected with >97.5% probability in B. 


coefficient of 2, 4, or 5/3 is applied, respectively. Of course, this re-normalization makes sense only for comparison, and should 

not be used otherwise. _ _ 

The results are shown in Fig. 10 and discussed in Sections VI-A to VI-C Each simulation is averaged on 10^ runs. All 


plots show AVAR (blue), MVAR (green) and PVAR (red) for the FP process, with the two-sigma uncertainty bars, and the SP 
process (grey). We set the reference value of the SP process at the lowest level that PVAR can detect with a probability of 
97.5%, i.e., at the upper point of the two-sigma uncertainty bar at r = 1024 s. This is highlighted by a grey circle at r = 1024 

s. 


A. Noise detection in the presence of white PM noise (Fig. \T^-b) 

White PM noise is a limiting factor in the detection of other noise processes because it is the dominant process in the 
front end of most instruments used to assess the frequency stability. We show the effect of white PM in two opposite cases, 
white FM noise and frequency drift. The former is present in all atomic standards, while the latter is present in all oscillators 
and standards, except in the primary standards. Frequency drift is a severe limitation in cavity stabilized lasers, and in other 
precision oscillators based on the mechanical properties of an artifact. 

The classical AVAR is clearly a poor option because of its response to white PM, versus the of the other variances. 
This is confirmed in our simulations. 

It is seen on Fig. p^A-B that in both cases MVAR cannot detect the slow process. The lowest value of MVAR (green plot) 
at 97.5% confidence (grey circle at r = 682 s) exceeds the reference grey line. 

The conclusion is that PVAR exhibits the highest detection sensitivity in the in the presence of white PM noise. 


B. Detection of flicker FM noise in the presence of white FM noise (Fig.\l^C) 

The detection of frequency dicker in the presence of white FM noise is a typical problem of passive atomic standards. Such 
standards show white FM noise originated from the signal to noise ratio, and in turn from beam intensity, optical contrast, 
or other parameters depending on the physics of the standard. Generally, after the white FM noise rolls off, cry(r) hits the 
dicker door. Cs clocks are a special case because they do not suffer from random walk and drift. So, dicker of frequency is 
the ultimate limitation to long-term stability, and in turn to timekeeping accuracy. In commercial standards, dicker FM exceeds 
the white FM at approximately 1 day integration time. Thus, fast detection of dicker enables early estimation of the long term 
behavior, and provides a useful diagnostic. 

We see on Figure [^C that the three variances show similar performances, with a small superiority of AVAR and PVAR. 
Again, MVAR suffers from the wider support of the wavelet, 3r instead of 2r. AVAR has a distinguished history of beeing 
the favorite tool of time keepers. 


C. Detection of slow phenomena (Fig. lOD-E-F) 


It is often useful to detect the corner where random walk or drift exceed the dicker door, or where the drift exceeds the 
random walk. This problem is typical of Rb clocks and H masers, and also of precision oscillators based on mechanical 
properties of a resonator. Our simulation shows that AVAR is superior, but PVAR has a detection capability close to AVAR. 
Conversely, MVAR is the poorest choice. 
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Fig. 10. Corner detection for the most common noise types (grey circles). The 97.5 % upper bounds of the confidence intervals over the variance estimates 
are figured by dashed lines (blue for AVAR, green for MVAR and red for PVAR). The lowest detected noise or drift by PVAR is represented as a solid black 
line. 


VIL Discussion and Conclusion 

PVAR is wavelet-like variance broadly similar to AVAR and MVAR, and intended for similar purposes. It derives from 
AVAR and MVAR after replacing the 11 and A counter with the Q counter, in turn based on the linear regression on phase 
data |[^. 

On closer examination, we notice that AVAR and MVAR address different problems. In the presence of white PM noise, 
MVAR has a dependence as 1/r^ instead of 1/r^. This is a good choice in microwave photonics and in other applications 
where the measurement of short term stability is important. The problem with MVAR is that the wavelet spans over 3r instead 
of 2r. Hence, AVAR is preferred for the measurement of long term stability and in timekeeping, where the largest value of r 
on the plot is severely limited by the length of the available data record. PVAR on the other hand is a candidate replacement 
for both because it features the 1/r^ dependence of MVAR and the 2r measurement time of AVAR. 

PVAR compares favorably to MVAR because it provides larger EDF, and in turn a smaller confidence interval. The objection 
that PVAR gives a larger response to the same noise level (right hand column of Table is irrelevant because the response is 
just a matter of normalization. It is only in the region of fast processes that AVAR has higher EDF than PVAR (Fig. [^, but 
this happens where AVAR is certainly not the preferred option. 

The best of PVAR is its power to detect and identify weak noise processes with the shortest data record. We have seen in 
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Sec. VI that PVAR is superior to MVAR in all cases, and also superior to AVAR for all short-term and medium-term processes, 
up to flicker FM included. AVAR is just a little better with random walk and drift. 

In conclusion, theory and simulation suggest that PVAR is an improved replacement for MVAR in all cases, provided the 
computing overhead can be accepted. Whether or not AVAR is preferable to PVAR for timekeeping is a matter of discussion. 
AVAR renders the largest r with a given data record. This is the case of random walk and drift. By contrast, PVAR is superior 
at detecting the frequency flicker floor, which is the critical parameter of the primary frequency standards used in timekeeping. 
These standards are supposed to be free from random walk and drift. Otherwise, when rendering the largest r is less critical, 
PVAR is until now the best option. 
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